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Abstract 

This paper develops and analyzes some continuous interior penalty finite element methods 
(CIP-FEMs) using piecewise linear polynomials for the Helmholtz equation with the first order 
absorbing boundary condition in two and three dimensions. The novelty of the proposed methods 
is to use complex penalty parameters with positive imaginary parts. It is proved that, if the 
penalty parameter is a pure imaginary number 17 with < 7 < C, then the proposed CIP-FEM 
is stable (hence well-posed) without any mesh constraint. Moreover the method satisfies the error 
(-^ I estimates Cikh + C2k'^h? in the H^-norm when k^h? < Co and Cikh + ^ when k'^h^ > Co and 

kh is bounded, where k is the wave number, h is the mesh size, and the C's are positive constants 
independent of k, h, and 7. Optimal order error estimates are also derived. The analysis 
is also applied if the penalty parameter is a complex number with positive imaginary part. By 
taking 7 0+, the above estimates are extended to the linear finite element method under the 
condition k^h^ < Co- Numerical results are provided to verify the theoretical findings. It is shown 
^ ■ that the penalty parameters may be tuned to greatly reduce the pollution errors. 

■ Key words. Helmholtz equation, large wave number, continuous interior penalty finite element 

I methods, prc-asyniptotic error estimates 
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^ . 

5 ■ 1 Introduction 

The problem of short waves (or waves with high wave numbers) in acoustics, electromagnetics or 
surface water wave applications was listed as an unsolved problem in finite element methods (FEMs) 
, in 2000 by Zienkiewicz in his review paper |38) . It still remains open although some big progresses 

' have been made since then. In this paper, we consider the following Helmholtz problem: 

(1) -Am - fc^M = / in rj, 

(2) — + iku^g on T, 

on 

where C R'', d = 2, 3 is a polygonal/polyhedral domain, F := dft, i = T denotes the imaginary 
unit, and n denotes the unit outward normal to dQ. The above Helmholtz problem is an approximation 
of the following acoustic scattering problem (with time dependence e'"*): 

(3) -Au ~k^u^ f in R'', 

(4) ^/ 9(7/ - u""=) _^ _ \ ^ ^ asr=|a;|^oo, 

\ or ) 

where m'"*^ is the incident wave and k is known as the wave number. The Robin boundary condition 
([2]) is known as the first order approximation of the radiation condition (j4]) (cf. I22]). We remark 
that the Helmholtz problem ([T])-([2]) also arises in applications as a consequence of frequency domain 
treatment of attenuated scalar waves (cf. 21 ). 
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The difficulties of FEMs applied to the Helmholtz problem ([I])-® with large wave number lie in 
both their theoretical analysis and numerical efficiency mainly due to the high indefiniteness of the 
problem. While for the one dimensional (1-D) case the FEMs have been well understood. Ihlenburg 
and Babuska [HI] proved that the linear FEM for a 1-D Helmholtz problem satisfies the following error 
estimate under the mesh constraint kh < 1. 

(5) ||V(w-ii^™)||^3(o) <Cikh + C2k'^h^. 

Here h is the mesh size and Ci,i = 1,2 are positive constants independent of k and h. Note that the 
first term on the right hand side of ([5]) is of the same order as the interpolation error in _ff ^-seminorm. 
It dominates the error bound only if k^h is small enough. The second term C2k{kh)^ dominates if 
kh is fixed and k is large enough. We remark that the condition of fixed kh, i.e., several points per 
wavelength, is sometimes used as the "rule of thumb" in the context of the numerical treatment of the 
Helmholtz equation. The estimate ([5]) says that this rule of thumb may give wrong results for large 
wave number k. The second term is called the pollution error of the finite element solution. In one 
dimension, the pollution effect can be eliminated completely by a suitable modification of the discrete 
bilinear form (cf. [6j^). However, the story for two and three dimensional Helmholtz problems is 
much different. It is shown that, in two dimensions, the pollution effect can be reduced substantially 
but cannot be avoided in principle (cf. [6l0]). As for the error estimates, to the best of the author's 
knowledge, no analysis for the linear FEM in two or three dimensions has been done when k^h is large. 
Note that the H^- and L^- error estimates can be derived by the so-called Schatz argument if k^h is 
small enough (cf. [9, 34 ) but this condition is too strict for large k. We refer to fJB] for a nice review 
on various FEMs for time-harmonic acoustics governed by the Helmholtz equation. For results on hp- 
FEMs, we refer to [32] but will not discuss here since we concern only methods using linear elements 
in this paper. The author would like to mention that, Engquist and Ying l24l [23] proposed recently 
some sweeping preconditioners for central difference schemes for the Helmholtz equation which have 
linear application cost and the preconditioned iterative solver (GMRES) converges in a number of 
iterations that is essentially independent of the number of unknowns or the frequency. Although the 
sweeping techniques are well possible (or have already been) applied to the linear FEM to provide 
efficient fast solvers, this combination is a not good candidate for efficient algorithm for the Helmholtz 
problem with large wave number, since the linear FEM itself is inefficient due to its (big) pollution 
effect. Next, we will not go any further on the issue of fast solvers and focus on the stability and error 
analyses of the schemes based on linear elements. 

In [25[ 126] ■ Feng and the author proposed and analyzed some interior penalty discontinuous 
Galerkin (IPDG) methods using piecewise linear polynomials for the problem ll])-© in two and 
three dimensions. It was proved that the proposed methods are unconditionally (with respect to mesh 
size h) stable and well-posed for all wave numbers A; > 0. Moreover, under suitable assumptions on 
the penalty parameters, the following error estimates were proved. 

\\u - Mhlli^ < Cikh + C2fc^/^/l*/^ if kh < 1, 

\\u - < C\kh + C2fc^/l^ if k^h^ < Co. 

where is some broken i7^-norm, A< B means A< C B, and the constants C's are positive and 

independent of fc, h, and the penalty parameters. Numerical tests show that it is possible to greatly 
reduce the pollution error and achieve better numerical results than FEMs by tuning the penalty 
parameters (see ^.25J). The discontinuous Galerkin (DG) methods (sometimes called discontinuous 
finite element methods) which are initiated in seventies of the last century (cf. [S] [71 H] [20] [37] [3] ) , use 
piecewise polynomials (or problem dependent functions) as trial and test functions. The continuity 
of the discrete solution across the interior edges/faces of elements is enforce weakly by introducing 
penalty terms or numerical fiuxes. As is well known now, DG methods have several advantages over 
the (continuous) FEMs, such as, local mass conservation, flexibilities in constructing trial and test 
spaces and meshes, additional parameters that may be tuned for some particular purposes. While one 
disadvantage is that a DG method usually has larger number of total degrees of freedom (DOFs) than 
the FEM. For example, on a given triangulation of il, the number of total DOFs of the linear IPDG 
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method is about six times of that of the hnear FEM in two dimensions and about 20 times in three 
dimensions. We refer the reader to [U [H [27l [28l [33l |39] and the references therein for other works on 
DG methods for Helmhohz problems. 

The purpose of this paper is to propose and analyze a linear continuous interior penalty finite 
element method (CIP-FEM) for the Helmhohz problem The CIP-FEM uses the same contin- 

uous piecewise linear finite element space as the linear FEM but modifies the sesquilinear of the FEM 
by adding a penalty term on the jumps of the flux across the interior edges/faces between elements, 
i.e., 



J(u, v) 



du 



dv 



where 7 > and £[ is the set of interior edges/faces. Note that the CIP-FEM in this paper uses a 
pure-imaginary penalty parameter i7 instead of a real one as the usual CIP-FEM does. This is helpful 
for theoretical analysis and numerical stability. It is should be remark that if the penalty parameter 
i7 is replaced by a complex number with positive imaginary part, the ideas of the paper still apply. 
Here we set its real part to be zero in the theoretical analysis for the ease of presentation. Let be 
the CIP finite element solution and let uf^^^ be the finite element solution. The following results are 
obtained. 

(i) The CIP-FEM attains a unique solution for any k > 0, h > and 7 > 0. 

(ii) There exists a constant Cq > independent of k, h, and 7, such that if fc > 1 and < 7 < 1, 
then the following stability and error estimates hold: 

( Cikh + C2k^h^, if k^h^ < Co, 



C^kh+—, if fc^/i^ > Co and yfc/i < 1, 

7 



where ^ (l|Vu||i2m) + 



(iii) Suppose k^h^ < Cq and fc > 1. Then the following estimates hold for the finite element solution 

„,FEM 

.FEM I 



|V(w-ur^)|U,„. < Cikh + C2k^h', 



(iv) Estimates in the L^-norm are also obtained. 

(v) Numerical tests show that the penalty parameters may be tuned to greatly reduce the pollution 
errors. 

The CIP-FEMs were originally proposed by Douglas and Dupont [50] for second order elliptic and 
parabolic problems and have been shown to have advantages for advection dominated problems |12[ 
[T51 im IT31 [IB] . Similar interior penalty procedures for FEMs utilizing continuous functions have also 
been introduced for biharmonic equations [71 1111 etc] . 

This paper is organized as follows. The CIP-FEM is introduced in Section [21 Some stability 
estimates are derived in Section [3l for any fc > 0, > 0, and 7 > 0. In Section |4l pre- asymptotic 
error estimates in H^- and L^-norms are proved for fc > 0, /i > 0, and 7 > by utilizing the error 
analysis for an elliptic projection, the stability results for the CIP-FEM, and the triangle inequality. 
In Section [5l the stability estimates in Section [3l and the error estimates in Section |4l are improved to 
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be of optimal order under the condition that k^h^ is small enough by using the technique of so-called 
"stability-error iterative improvement" developed in [26 . In Section [6l the well-posedness, stability 
and error estimates for the linear FEM are established under the condition that k^h^ is small enough 
by taking the limits of the estimates for the CIP-FEM as the parameter 7 — > 0-1- . 

Throughout the paper, C is used to denote a generic positive constant which is independent of 
h, k, and the penalty parameters. We also use the shorthand notation A < B and B > A for the 
inequality A < CB and B > CA. ^ ~ S is a shorthand notation for the statement A < B and 
B < A. We assume that k > 1 since we are considering high-frequency problems. For the ease of 
presentation, we assume that k is constant on the domain fl. 



2 Formulation of continuous interior penalty finite element 
methods 

To formulate our CIP-FEMs, we first introduce some notation. The standard space, norm and inner 
product notation are adopted. Their definitions can be found in [inillH]- In particular, (•, •)q and (•, 
for E C dQ denote the L^-inner product on complex- valued L^{Q) and iv^(E) spaces, respectively. 
Denote by (•, •) := (•, On and (•, •) := (•, ■)g^. 

Let 7^ be a family of triangulations of the domain fl parameterized by /i > 0. For any trian- 
gle/tetrahedron K G Th, we define Hk ■— diam(ii'). Similarly, for each edge/face e of if € Th, define 
he := diam(e). Let h = maxKeTh ^k- We assume that the elements of Th are shape regular. We 
define 

£[ := set of all interior edges/faces of Th, 

:~ set of all boundary edges/faces of Th on F. 

We also define the jump [v] of v on an interior edge/face e = dK n dK' as 

r , I J v\k — v\k', if the global label of K is bigger, 

\ v\k' — v\k , if the global label of K' is bigger. 

For every e = dK D dK' g f ^, let Ue be the unit outward normal to edge/face e of the element K if 
the global label of K is bigger and of the element K' if the other way around. For every e € £^ , let 
rie = n the unit outward normal to 

Now we define the "energy" space E and the sesquilincar form a/i(-, ■) on E y. E as, follows: 

E■.= H^{n)c^ W H^{K), 

KeTh 

(6) a/i(u, 11) := (Vm, Vf) -I- J(u, w) \/u,v^E, 



where 
(7) 



J{u, v) 



E 



du 



dv 



and "fejS & £h are nonnegative numbers to be specified later. 

Remark 2.1. (a) The terms in J{u, v) are so-called penalty terms. The penalty parameter in J{u, v) 
is ije ■ So it is a pure imaginary number with positive imaginary part. It turns out that if it is replaced 
by a complex number with positive imaginary part, the ideas of the paper still apply. Here we set their 
real parts to be zero partly because the terms from real parts do not help much ( and do not cause any 
problem either) in our theoretical analysis and partly for the ease of presentation. On the other hand, 
our numerical experiments in Section^ indicate that using penalty parameters with nonzero real parts 
helps to reduce the pollution effect in the error. 
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(b) Penalizing the jumps of normal derivatives was used early by Douglas and Dupont 120] for 
second order PDEs and by Babuska and Zldmal /?/ for fourth order PDEs in the context of finite 
element methods, by Baker ^8] for fourth order PDEs and by Arnold I3j for second order parabolic 
PDEs in the context of IPDG methods. 

(c) In this paper we consider the scattering problem with time dependence e'"*, that is, the signs 
before i '5 in the Sommerfeld radiation condition Q and its first order approximation ^ are positive. 
If we consider the scattering problem with time dependence e~"^*, that is, the signs before i's in (|4]) 
and ([2]) are negative, then the penalty parameters should be complex numbers with negative imaginary 
parts. 



It is clear that J{u, v) 
II])-©, then 



if M e H'^lfl) and v e E. Therefore, if u e H'^{Q) is the sohition of 



(8) aft(u,w) - fc^(u,u) + ifc(M,u) = + , Vv e E. 
Let Vh be the hnear finite element space, that is, 

Vh := {vh e H\n) : VhU e Pi{K), € %} • 

where Pi{K) denote the set of all linear polynomials on K. Then our CIP-FEMs are defined as follows 
: Find u/i S Vh such that 

(9) ah{uh,Vh) - k^{uh,Vh) + ik{uh,Vh) = {f,Vh) + {g,Vh) , Vv/i € V?,. 
The following semi-norm on the space E is useful for the subsequent analysis: 



(10) 



\l.h 



\\Iv\ 



L2(0) 



dv 



dUe 



1/2 



In the next three sections, we shall consider the stability and error analysis for the above CIP- 
FEMs. Especially, we are interested in knowing how the stability constants and error constants depend 
on the wave number k (and mesh size h, of course) and what are the "optimal" relationship between 
mesh size h and the wave number k. For the ease of presentation, we assume that 7e ~ 7 for some 
positive constant 7 and that hx — h. 



3 Stability estimates 

We first recall the stability estimates for the original Helmholtz problem ([T])-([2|) (cf. [I9l 130) ). 

Theorem 3.1. Suppose C R'' is a strictly star-shaped domain. Then the solution u to the problem 
satisfies 

(11) MHHn)<k'-'{\\f\\mn) + \\9\\mr)) 

for j ^0,1 ifue i?3/2+£(f]) jor some £ > 0. (HH) also holds for j = 2 if u E H'^{n). 

The key idea in their analysis is to test 1^ hy v — u and v = {x — xq) • Vw, respectively, and 
use the Rellich identity (for the Laplacian), where xn is a point such that the domain ft is strictly 
star-shaped with respect to it. The idea has been successfully applied to the discontinuous Galerkin 
methods (cf. [25l [H |27j ) and to the spectral-Galerkin methods (cf. [3F). As for our CIP-FEMs 
although the test function Vh = uu can still be used, the test function Vh = {x — xq) ■ \/uh does not 
apply since it is discontinuous and hence not in the test space Vh . For stability results for other types 
of boundary conditions we refer to [ITl [32] . 

Next, we derive stability estimates for the CIP-FEMs Note that Uh is piecewise linear on 
Th and hence Auh — on each element K G 7;, . We will show, by using integration by parts 
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elementwisely, that || Vu;i||^2(f^-) may be bounded by the jumps of across each interior edge/face 
e e and the L'^{n) -norm and the L^(r)-norm of u^. Moreover the coefficient before ^2(5-2-) 
can be controUed. On the other hand, by taking the test function Vh — Uh in ([9]), we may derive 
some reverse inequahties, that is, bound the jumps of across e G and the norms of Uh by 



L2(n) 



and the given data. Then the desire stability estimates follow by combining them. 



We first bound ||Vu/i||^2(s^-) by using integration by parts on each element. 
Lemma 3.2. For any < e < 1, there exists a constant such that 



Proof. Noting that Uh is piecewise linear, we have 

|2 



^(r) ^ p/j2 

' edS' 



dUh 
dUe 



IVuJ 



L2(n) 



Ken 



KeTh 



duh 
dK dn 

duh 



Uh 



dUe 



Uh 



< E W^^hWme) hh\\L2(e) + E 



duh 

dUe 



l"/i|lL2(e) 



For any edge/face e € £h, let Kg e 7/i be an element containing e. From the trace inequality and the 
inverse inequality. 



I Vu,i||^2(o) <C Y h^ iiy^^ii ^^^^^^ 11^^11 ^^^^^ +(7^/1/ 



/2 



eeSi 



< Ch-'/^ WVuhh^^n) IWhh^^r) + C^-^'^h-^ ( 7e/ie 

C 

< £ ||V'«/i||i2(jj) + ||uh|li2(r) 



duh 

dUe 



duh 

2 

LHe) 



1/2 



"'illL2(_ftr,) 



"'i|lL2(0) 



+ £(l-e)/c2||u^||2 , 



C 



^ ' ee£' 



E ^'^'^ 



LHe) 



which implies that Lemma 13.21 holds. 

Then we derive some reverse inequalities by taking Vh = u/i in ([9]). 
Lemma 3.3. Let Uh G Vh solve ([9]). Then, 



□ 



(12) 
(13) 



E '^^'^'^ 



duh 



k \\uh 



fc2 
2 

lL2(r) 



C 



C 



Proof. Taking Vh ~ Uh in ^ yields 



(14) 



ah{uh, Uh) - fc^ ||uh|li2(n) + ifc Il'«/i|li2(r) = (/. "/i) + (5, ""/i) • 



L2(r) 
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Therefore, by taking real part and imaginary part of the above equation we get 



(15) 
(16) 



hh\\l2(^n) - l|Vuh||L2(o) < \if,uh) + {g,uh)\, 



duh 



drie 



+ k ||uh|li2(r) < \{f,Uh) + {9,Uh)\ 



From (H 



duh 



dUe 



1 k 

+ k ||uh|li2(r) < II/IIl^(o) hftllL2(J7) + ^ hWl^D + 2 ll"''lli^(r) 



which imphes 

(17) ^-f"- 



duh 



k 1 

L2(e) ^ ^'^ 



On the other hand, from (|T5]) . 



1 fc 

l|M/.||i2(o) < ||Vm,,||^2(o) + WWL^iu) + ^ Il5lli2(r) + ^ Il"'>lli2(r) 

By combining the above two estimates, we conclude that 



||uh|li2(iy < ||Vuft||^2(o) + 2 ||/||£2(o) + T Il3lli2(r) 



2 fc^ 



ll3llL2(r) 



which implies that ([T2|) holds. 

Plugging (|12p into the right hand side of ([TTl) yields 



L2(e) 



, ^ II l|2 
+ ^ l|t^/i|lL2(r) 

1 



I|Vu/i|Il2(o) + ^ ll/llL^(i2) + \\9\\l^(T) 



^llffllW) 



c c c 

< -j: WfWh^n) l|Vwh||L2(o) + p 11/11^2(0) + Il5lli2(r) 
That is, ([13]) holds. This completes the proof of the lemma. 



□ 

By combining Lemma 13.21 and Lemma 13.31 we may derive the following stability estimates for the 
CIP-FEMs. 

Theorem 3.4. The solution Uh G Vh to the scheme ((9]) satisfies the following stability estimates. 



(18) 
(19) 



ll^"'»llL2(r2) + k^ h'illL2(f^) < C^t^ \\f\\\--{n) + C'sta Ilfflli2(r) 



duh 

dUe 



LHe) 



U II l|2 - ^' 

W^hW L2(r) 



<^\\f\\Un,^-,h\ 



2 

LHT) ■ 



Here 



C'sta • — 



1 



1 
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Proof. By taking £ = - in Lemma 13.21 and applying Lemma 13.31 



\yuh\ 



2 ^ ; 2 II ||2 

<^k \\uh 



^^(") ^ kh 



C 

— k\\uh\\\2(j'^ 



C 



< 



^(2||Vu^ 



1l2(J2) 



p WfWL^n) 



C II ||2 

J \\9\\LHr) 



duh 



kh k^h'^"f 



1 



l/ll 



1 



^2 ll/llL2(f2) 



2 

lL2(r) 



2„ „2 
<3 l|Vw„"' 



1 

' kh " 



CI 



-(- 

k \kh 
1 



k^h"^"/ 



/llL2(r!) ll'^'"''liL2(n) 



k^h^-fJ Vfc2 



l/ll 



1 



L2(ri) 
1 



1 II l|2 

rll-^lli^cr) 



1 



k'^h^^J k 



fc2ft,2^ 
rll5lli2(r). 



2 1 



II/IIl2(o) 



Therefore, 



II^"'i|Il2(J1) ^ C's^ta II/IIl2(o) + C'sta ||.9llL2(r) ■ 

Then the proof of the theorem foUows by combining the above estimate and Lemma [ 
Corollary 3.5. The CIP-FEM (|9]) attains a unique solution for any k > 0, h > and 7 > 0. 



□ 



Remark 3.1. (a) For the general case when the penalty parameters or the meshes may be nonuniform, 
Theorem \3.4\ and Corollaru \3. 5\ still hold with 7 and h replaced byj = min^ggj andh = minxeTh f^K, 
respectively. The proof is similar and is omitted. 

(b) If J ^ 1 ^u9 ihen Csta 1 which implies the following stability estimates for the CIP-FEM: 



k^h^ 



1 



IWhWi^h ^ II/IIl2(o) + ll5llL2(r) ' IWhh^in) ^ liWfh^n) + ll.9llL2(r) )■ 

These estimates are of the same order as those for the Helmholtz problem ((T])-([5]) (cf. Theorem \3.1]) . 
But we do not suggest to choose 7 as above when k'^h^ is small, since a large 7 may cause a large 
error of the discrete solution (cf. Theorem \4.4\ below). 

(c) The stability estimates in Theorem \S.4\ will be improved in Section\^when k'^h^ is small enough. 

Note that if k'^h^ > 1 and 7 < 1, then Csta ^ — arid hence 

7 

IWf^hn < ^{Wfhnn) + Mlhd), IWkh.^n) < ^( II/IIl^(o) + ll.9llL^(r) ) ■ 



4 Error estimates 

In this section we first introduce an elhptic projection of the solution u to the Helmholtz problem ([1])- 
^ and estimate the error between them. Then we estimate the error between the elliptic projection 
and the CIP finite element solution by using the stability estimates in the previous section. In 
what follows, we assume that the domain f2 is a convex polygon/polyhedron. Then u G II^{Q.) (cf. 
and Theorem 13.11 implies that 

(20) ||ii||^.(,,)</cM(/,5), 



where M{f,g) = \\f\\L2(a) + WaWmr) 
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4.1 Elliptic projection and its error estimates 

For any w ^ E, we define its elliptic projection Wh € Vh by 

(21) ah{wh,Vh) + ik{wh,Vh) ^ ah{w,Vh) -\-ik{w,Vh) Vw/i G V^. 

In other words, Wh is an CIP finite element approximation to the solution w of the following (complex- 
valued) Poisson problem: 

-Aw = F in f], 

dw 

— (- ikw = ^ on r, 

on 

for some given functions F and which are determined by w. 

The following lemma gives the continuity and coercivity of the sesquilinear form ah{-, ■) whose 
proof is obvious and is omitted. 

Lemma 4.1. For any v,w ^ E, 

(22) \ah{v,w)\, \ahiw,v)\ <\\v\\^jj\w\\^j^, 

(23) Rcah{v,v) + lmah{v,v) \\v\\l^^ . 

Let u be the solution of problem (HJ-© and Uh be its elliptic projection defined as above. Then 
(j2ip immediately implies the following Galerkin orthogonality: 

(24) ah{u-Uh,Vh) +ik(u-Uh,Vh) =0 Vw;, € V^. 
Lemma 4.2. There hold the following estimates: 

(25) \\u-u^\\,,^<{l + j + khy^hhM{f,g), 

(26) \\u-uhh2^n) ^ {l + l + kh)kh^M{f,g). 

Proof. Let iih g Vh be the Pi-conforming finite element interpolant of u on the mesh Th- Then Uh 
satisfies the following estimates (cf. [TPl 118) ): 

(27) h^^h\\L2(_^) <h^\u\^^^^y l|V(u-Uft)||i2(f2) < 



which imply that 

(28) lk-Uft|lL2(r) < 

(29) ^ (l + 7)'^'/ih|ff2(o): 

where we have used the trace inequality ||w||^2(r) II''^'IIl2(I2) derive (|28p and used the 

local trace inequality ||w||^2(gjc) ^ ^k^"^ \\''^\\l'^(k) ^ ^]<^ \\'^''^\\l'^{k) ^'^^ K € Th io derive ([29)) . 
Let rj -.^ u — Uh- From (P^ . 

(30) a/i(ry, ?7) + ifc (77, 77) = ah{r], u - Uh) + ik (77, u- Uh) ■ 
It follows from Lemma 14.11 and (l30l) that 



kllL =Reah(ry,?7) +Ima,i(?7,?7) 
= Re(ah(?], ?7) + ifc (?7, ry)) 

+ Im (a/i(?7, 77) + ifc (77, rj)) - k (77, ry) 

= Re(a/i(7], M- Mh) +ifc wh)) - fc||f?|lL2(r) 

+ Im(a/i(77,u - Uh) + ifc (t?, u - Uh)) 

<c{\\r]\\^j, \\u - iihWi.h + ^ ll^llL2(r) h - "/^^^(r)) " ^ ll'^lli^cn 
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Therefore, it follows from (gH), (HH), and (gHl) that 

<{l + -/ + kh)k^h^M{f,gf 



(31) hllL + hllL2(r) ^ h - UhWl^h + k\\u- Uh|li2(r) 



That is, ^ holds. 

To show (|26l) . we use the Nitsche's duality argument (cf. [10l[T8]). Consider the following auxiliary 
problem: 

(32) —Aw — u — Uh in ri, 

duo 

— ikw — Q on r. 

on 

It can be shown that w satisfies 

Let Wh be the Pi-conforming finite element interpolant of w on 7h. Testing the conjugated of (j32]) by 
u — Uh and using ([M)) we get 

\\u - '«/i|lL2(f^) = -(u - Uh, Aw) = ah{u - Uh, w) + ik{u- Uh,w) 
= ah{u - Uh, w - Wh) + ik{u- uh, w - lih) 
^ II" - UhWiM 11^ - ^hWi^h + k\\u- M/i||i2(r) \\w - Wh\\L2(^r) 

which together with ([31]) and (1551) gives (PB)) . The proof is completed. □ 
4.2 Error estimates for the CIP-FEMs 

In this subsection we shall derive error estimates for the scheme ([9]). This will be done by exploiting 
the linearity of the Helmholtz equation and making use of the stability estimates derived in Theorem 
13.41 and the projection error estimates established in Lemma [4.21 

Let u and Uh denote the solution of ([l])-([2]) and that of ([9]), respectively. Define the error function 
eh := u — Uh- Subtracting ([9]) from ([8]) with v = Vh & Vh yields the following error equation: 

(34) ah[eh,Vh) -k^ieh^Vh) + ik{eh,Vh) ^Vh^Vh- 

Let Uh be the elliptic projection of u as defined in the previous subsection. Write eh — V ~ C with 

ry.^u-Uh, £,:=Uh-Uh- 

From and ([24]) we get 

(35) tthit Vh) - Vh) + ik (^, Vh) = 0^(77, Vh) - k^{ri, Vh) + ik {rj, Vh) 

= -k^{ri,Vh) \fvh e Vh. 

The above equation implies that ^ e V/i is the solution of the scheme (O with source terms / — — fc^ry 
and g = 0. Then an application of Theorem 13.41 and Lemma 14.21 immediately gives the following 
lemma. 

Lemma 4.3. ^ ^ Uh ~ Uh satisfies the following estimate: 

(36) llelli,, + k Uh-in) ^ Cst.{l + 7 + kh)k^hHl{f, g), 
where Csta is defined in Theoreni\3.4\ 
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We are ready to state our error estimate results for scheme ©, which follows from Lemma [4.21 
Lemma 14.31 and an application of the triangle inequality. 

Theorem 4.4. Let u and Uh denote the solutions of and respectively. Then there exist 

two positive constants Ci and C2 such that the following error estimates hold. 

(37) \\u-uh\\^j^< {l+j + kh){Cikh + C2Cst.k''h^)M{f,g), 

(38) \\u-Uh\\L2(^n) < {l+l + kh){Cikh^ + C2Cst.k^h^)M{f,g), 

where Cgta is defined in Theorem \3.4\ and AI{f,g) — ||/||i2(-j^') + llffHiaj-p-) • 

Remark 4.1. (a) If kh < 1, k^h^ > 1, and 7 < 1, then Csta ^ ,0,0 and we have the following error 
estimates for the CIP-FEM: 

II" - M\i,h ^ {Cikh + ^)Mif,g), \\u ~ u4^.^^^ < {Cikh^ + ^)M{f,g). 

The pollution term in the above error estimate is 0(1) if j ~ 1. By contrast the pollution term for 
the linear FEM when k^h^ > 1 is expected to be of order k^h^ as that proved for the one dimensional 
case (cf. JSTJ). 

(b) The error estimates in Theorem \4-.4\ will be improved in the next section when k^h^ < Cq for 
some constant Cq independent of k, h, and the penalty parameters. 

5 Stability-error iterative improvement 

In this section we improve the stability estimates in Theorem 13.41 and the error estimates in Theo- 
rem |4j4] under the condition that k^h^ is small enough, by using the trick of so called "stability-error 
iterative improvement" developed in [26]. 

Theorem 5.1. Letu anduh denote the solutions of and respectively. Assume that"/ < 1. 

Then there exists a constant Cq > 0, which is independent of k, h, and the penalty parameters, such 
that if k^h^ < Cq, then the following stability and error estimates hold: 



(39) h/Jli,„<M(/,ff), 

(40) \\nh\\mn)<lMif,9), 

(41) \\u~Uh\\^,^<{Cikh + C2k^h^)M{f,g), 

(42) h~"/>.|lL.(o) ^ {Cikh' + C2k'h')M{f,g), 



where M{f,g) = ||/||i2(o) + ll5llL2(r)- 

Proof. It suffices to prove ([5^ . since PO]) follows then from Lemma [3.31 (specificallv. (UH)) and (jH))- 
follow from the improved stability estimates and the argument used in the proof of Theorem l4.4l 
Suppose Csta > 1, otherwise, ([M)) holds already (cf. Theorem 13. 4p . 
From Theorem 13.41 we have, for any / S L^(r2) and g S L'^{r), 

(43) ll"^.|li,,.<C^staA/(/,5). 
Suppose kh <1. Then (|35|) and Lemma 221 imply that 

ll"'i - "'illiji ^ C'stafc^ II"- "/i|Il2(J7) ^ Cst!,k^h^M{f,g). 
Therefore from the triangle inequality and Lemma 14.21 we have 

(44) \\u-Uh\\i,^ < \\u-Uh\\ij, + \\uh-Uh\\ij^ < {kh + Cst^k^h'^)Mif,g). 
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Now it follows from the triangle inequality and Theorem 13.11 that 

(45) WuhWi^h < ll^lliji + \Wh-u\\^^^ = ||Vu||^2(o) + 

\ + kh + at.k''h^y4{f,g). 



< 



Repeating the above process yields that there exists a constant C independent of k, h, and the penalty 
parameters, and a sequence of positive numbers Aj such that 

(46) \\u„\\,,^<A,M{f,g), 
with 

Ao-Csta, A, =Cil + kh)+Ck^h^A,^i, j = l,2,---. 

A simple calculation yields that ii C k^h^ < 6 for some positive constant 9 < 1 then 

, C{l + kh) 
hm A, = ^ -, 

which implies (|39|) . □ 

Remark 5.1. (a) Note that the stability estimates in (|39p and (j40p are of the same order as the PDE 
stability estimates given in Theorem \3.1\ 

(b ) Note that the estimates in Theorem \5.1\ are uniform with respect < 7 < 1 . In the next 
section, by passing to the limit 7 — > 0+ in the CIP-FEMs ([9|) and in (|39p - (I42p , we will give stability 
and error estimates for the FEMs. 

6 Stability and error estimates for the linear finite element 
method 

It is clear that both the bilinear form a/i(-, •) and the CIP finite element solution Uh to ([9]) depend 
on the penalty parameters 7e. In this section, we choose 7e = 7 and denote by a^(-, •) :— ah{-, •) and 
by uj^ := Ufi. Obviously, if 7 vanishes, the CIP-FEM ^ "degenerates" to the standard linear FEM: 
Find it™M g y^^ g^(,j^ ^Y^^^ 

(47) iVuf"^, Vvh) - fc2(«FEM^ ^ ^^FEM^ ^ ^ f^^^ ^ y^^ 

Next we will derive stability and error estimates for the FEM by showing that converges as 

7 0+. 

Theorem 6.1. There exists a constant Cq > independent of k and h such that if k^h^ < Co, then 
([T7|) attains a unique solution uj^^^ E Vh which satisfies the following stability and error estimates: 

(48) \\yur''\\r^.^n)^Mif,g), 

(49) |hr^lL.(^)<iAf(/,.9), 

(50) \\V{u-ul'''^)\\^ < {Cikh + C2k'h^)M{f,g), 



Il2(o) 

where M{f,g) = 11/11^2(0) + Uh^r)- 



(51) ll"-<™IL2(o) ^ {Cikh^ + C2k^h^)AI{f,g), 
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Proof. Let Cq be the constant defined in Theorem 15. II Suppose k^h^ < Cq and h is fixed. Note that 
the space Vh is finite dimensional and hence any two norms on Vh are equivalent. It is clear that, 
if uj^ converges to some function in H^{n) as 7 — ^ 0+, then the existence of a finite element 

solution to (H71) and the estimates for the finite element solution follow by letting 7 0+ 

in the CIP-FEM © and in Theorem 15. II Next we prove the convergence of u1 by using the Cauchy's 
convergence test. 

By letting 7 = 71, 72 in ([9]), respectively, and taking the difference, we get 

al {ul , z;,) - {uf , vh) - k\ul - uI\vh) + ik {u^ - , vh) = 0, e F,. 

Recall that 

aliv, v) = (V(^, Vv) + i7 ^ /le 

Clearly, uj^ — uj^ is the solution of the following discrete problem: 

- k\ul - ul\vu) + \k {ul - ul\vu) 
= (71 - 72)(A, Vh) := -i(7i - 72) X! 

Therefore, from Theorem 15. II 

(52) h7;-<1Ui(o)^|7i-72|||ML.(o)- 
Since any two norms on Vh are equivalent, from Theorem 15. II we have 

WUl^^) < C{h) \\uf\\H^(n) ^ C{h)M{f,g), 

where C{h) is some constant which is dependent of h but independent of 71 and 72. By combining 
the above two estimates, we have 

\WT~<lmm^b^-^^\C{h)M{f,g). 

Thus uj^ converges in H^{rt) as 7 ^ 0+. 

It remains to prove the uniqueness. Let be any solution to the FEM By repeating the 

lines for deriving (j52l) with 72 = 0, we obtain, 

hi' -^h\\m{n) ;$7i|IMIl2(o), 

where fh depends on but is independent of 71. Therefore lim^j_).ow7^ ~ u^-^ which implies that 
= '"/i^'^- This completes the proof of the theorem. □ 



dip 




dv 




5 


drie 



\fLp, V £ E. 







dvh 


dUe 


7 





yvh e Vh- 



7 Nmerical examples 

Throughout this section, we consider the following two-dimensional Helmholtz problem: 

(53) -Au-k^u^f:^^^^^^ inn, 

r 

Ou 

(54) — +iku^g on Tr := d^l. 

on 

Here ft is the unit regular hexagon with center (0, 0) (cf. Figure [1]) and g is so chosen that the exact 
solution is 

, cos(fcr) cos fc + i sin /c , , 

^"'^ " = — — kimrum'^^^ 
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Figure 1: Geometry (left) and a sample mesh Ti/j that consists of congruent and equilateral triangles 
of size h = 1/7 (right) for the example. 

in polar coordinates, where Jii{z) are Bessel functions of the first kind. 

Wc remark that this problem has been computed in [25] by using interior penalty discontinuous 
Galerkin methods. We use the same example for the convenience of comparison. 

For any positive integer m, let 7i/m denote the regular triangulation that consists of 6to^ congruent 
and equilateral triangles of size h = 1/m. See Figure [T] (right) for a sample triangulation Ti/y. We 
remark that the number of total DOFs of the CIP-FEM on the triangulation Ti/m is 3to^ + 3to + 1 
which is the same as that of the linear FEM and about one sixth of that of the linear IPDG method 
(cf. [25]). 



7.1 Stability 

Given a triangulation Tin recall that Uh denotes the GIF finite element solution and uj^^ denotes 
the Pi-conforming finite element approximation of the problem (|53 p - (|54p . In this subsection, we use 
the following penalty parameters for the CIP-FEM (cf. ([9])): 



(56) 



le 



= 7 = 0.1 Veef,^. 



Then, according to Theorem 13.41 and Theorem 15.11 we have the following stability estimate for the 
GIF finite element solution Uh- 



(57) 



where C„ta = t- 



L2(0) " 
1 



"'i|lL2(n) 

10 



< 



™ { + ll5llL2(r) , Csta ll/llL2(n) + ll5llL2(r) } 

Noting that the stability estimate in L^-norm is a direct consequence 



k ' k^h ' k^K^ 
of that in _ff^-seminorm (cf. Lemma 



we only examine the stability estimate for HVw/iH^s 



FigureOplots the iJ^-seminorm of the GIF finite element solution || Vu/ij 



the _ff ^-seminorm of 



the finite element solution || Vm™^||^2(q-) for h = 0.005 and 0.002, respectively, and the _ff^-seminorm 



of the exact solution || Vu||^2(f2) , for k 



1, • • • , 500. It shows that II Vm||^2(o) — 1, II Vw/i||i2(Q) 



< 1 and 

II Vuh|lL2(Q) decreases for k large enough as indicated by ([57|. It is also shown that || Vu™^||^2(f^) ^ 1 

which is guaranteed theoretically only for k^h^ small enough (cf. Theroeni l6.1|) . 

Figure O shows stability behaviors of Uh and when kh — 1 for k = 1,2,--- ,500, which 

satisfies the "rule of thumb". Note that Csta = 12/fc for kh = 1. It is shown that || Vu;j||^2(q) is in 
inverse proportion to k for k large which means the term Csta II/IIlsj-q-) dominates the stability bound 
for k < 500 (cf. (|57)) ). As a matter of fact, numerical integrations show that ||/||i2(Q) is more than 
25 times ||g||i2(r) for k = 200, • • • , 500. 
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Figure 2: ||Vu,i||^2(o) (solid), || Vu)^™||^2(^j (dashed) for h = 0.05 and 0.005, respectively. The 
dotted line gives the ff^-seminorm of the exact solution ||Vm||^2(q)- 
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k 

Figure 3: ||Vu^||^2(q) (solid) and c/k (dash dotted) versus to k with kh = 1, where c is some constant. 
The dashed line gives ||Vw)j^'^||^2(f2)- The dotted line gives the i/^-seminorm of the exact solution 
l|Vu||i2(o)- 

7.2 Error estimates 

In this subsection, we use the same penalty parameter (7 = 0.1) for the CIP-FEM as in ([56| . From 
Theorem 14.41 (cf. Remark 14.11) and Theorem 15.11 the error of the CIP finite element solution in the 
iJ-'^-seminorm is bounded by 

(58) |lV(M-Uh)|lL2(o) <Cikh + C2inm{k^h^,l} 

for some constants Ci and C2 kh < 1. On the other hand, from Theorem 16. 1[ the error of the finite 
element solution in the if ^-seminorm is bounded by 

(59) \\y{u - <™)||i2(o) < ^1^^ + ^2fc^/i' 

for some constants Ci and C2 if k^h^ < Cq. The second terms on the right hand sides of ([58|) and ([59|) 
are the so-called pollution errors. We now present numerical results to verify the above error bounds. 

In the left graph of Figure 21 the relative error of the CIP finite element solution with parameters 
given by (j56p and the relative error of the finite element interpolant are displayed in one plot. When 
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the mesh size is decreasing, the relative error of the CIP finite element solution stays around 100% 
before it is less than 100%, then decays slowly on a range increasing with fc, and then decays at a 
rate greater than —1 in the log- log scale but converges as fast as the finite element interpolant (with 
slope —1) for small h. The relative error grows with k along line kh = 0.25. By contrast, as shown 
in the right of Figure HI the relative error of the finite element solution first oscillates around 100%, 
then decays at a rate greater than —1 in the log- log scale but converges as fast as the finite element 
interpolant (with slope —1) for small h. The relative error of the finite element solution also grows 
with k along line kh = 0.25. 



10" 



10"' 











k = ri(p- 











10" 



10"' 























k = V 



l/ll 



l/h 



Figure 4: Left graph: the relative error of the CIP finite element solution with parameters given 
by (|56p (solid) and the relative error of the finite element interpolant (dotted) in Tf^-seminorm for 
fc = 5, fc = 10, fc = 50, and k — 100, respectively. Right graph: corresponding plots for finite element 
solutions. 

Unlike the error of the finite element interpolant, both the error of the CIP finite element solution 
and that of the finite element solution are not controlled by the magnitude of kh as indicated also 
by the two graphs in Figure [3 It is shown that when h is determined according to the "rule of 
thumb" , the relative error of the CIP finite element solution keeps less than 100% which means the 
CIP finite element solution has some accuracy even for large k, while the finite element solution is 
totally unusable for large k. Figure [5] displays the surface plots of the real parts of the linear interpolant 
of the exact solution (left), the CIP finite element solution with parameters given by ([56]) (center) , 
and the finite element solution (right), for k = 100 on the mesh with mesh size h = 1/100. It is shown 
that the CIP finite element solution has a correct shape although its amplitude is not very accurate. 
By contrast, the finite element solution has both wrong shape and amplitude. We remark that the 
accuracy of the CIP finite solution can be further greatly improved by tuning the penalty parameter 
i7, see Subsection 17.31 below. 

Next we verify more precisely the pollution terms in ()58p and ()59|) . To do so, we introduce the 
definition of the critical mesh size with respect to a given relative tolerance. 

Definition 7.1. Given a relative tolerance e and a wave number k, the critical mesh size h{k,e) with 
respect to the relative tolerance e is defined by the maximum mesh size such that the relative error of 
the CIP finite element solution (or the finite element solution) in H^-seminorm is less than or equal 
to e. 

It is clear that, if the pollution terms in ([58)1 and (l59l) are of order k^h^, then h{k,e) should be 
proportional to k~^/'^ for k large enough. This is verified by Figure [7] which plots h{k,0.5) versus 
k for the CIP finite element solution (left) with parameters given by (|56p and for the finite element 
solution (right), respectively. We remark that the maximum wave number such that h{k, 0.5) > 0.001 
is fcmax = 266 for the CIP-FEM with parameters given by ([5S)) and is fcmax — 280 for the FEM. 
Note that if the mesh size h — 0.001, then the number of total DOFs of the CIP finite clement 
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k k 

Figure 5: The relative error of the CIP finite eleraent solution (left) with parameters given by 
and that of the finite element solution (right) in if ^-seminorm computed for k — 1,2, ■ ■ ■ , 500 with 
mesh size h determined hy kh = 1 and kh — 0.5, respectively. 



system is 3, 003, 001, so is that of the FEM. Therefore, k — fcmax is the maximum wave number such 
that the problem ((55)) - ([5^ can be approximated by the CIP-FEM (or FEM) with relative error in 
iJ^-seminorm < 50% while using at most 3,003,001 total DOFs. 

7.3 Reduction of the pollution effect 

In |25J, it is shown that appropriate choice of the penalty parameters can significantly reduce the 
pollution error of the symmetric IPDG method. In this subsection, we shall show that the same thing 
holds true for the CIP-FEM. We use the following parameters: 

(60) i7e = 17 = -0.07 + O.Oli. 

We remark that this choice of i7e is the same as the choice of the penalty parameters i7i_e from [251 
Subsection 6.4] for the IPDG method. 

The relative error of the CIP finite element solution with parameters given by ([50]) and the relative 
error of the finite element interpolant are displayed in the left graph of Figure HI The CIP-FEM with 
parameters given by (IBH]) is much better than both the CIP-FEM using parameters given by ([55)) 
and the FEM (cf. Figure 2] and Figure [S]) . The relative error does not increase significantly with the 
change of k along line kh = 0.25 for k < 100. But this does not mean that the pollution error has 
been eliminated. For more detailed observation, the relative errors of the CIP finite element solution 
with parameters given by (|60p . computed for all integer k from 1 to 500 for kh = 1 and kh = 0.5, are 
plotted in the right graph of Figure [H It is shown that the pollution error is reduced significantly. 

Figure [S] plots h{k, 0.5), the critical mesh size with respect to the relative tolerance 50%, versus k 
for the CIP-FEM with parameters given by (pO]) . We recall that h{k,0.5) is the maximum mesh size 
such that the relative error of the CIP finite element solution in ff ^-seminorm is less than or equal to 
50%. The decreasing rate of h{k, 0.5) in the log-log scale is less than —1.5 for k from 1 to a relatively 
large value, which means that the pollution effect is reduced. We remark that the maximum wave 
number under the condition h{k,0.5) > 0.001 is fcmax ~ 622 for the CIP-FEM with parameters given 
by (|60l) which is more than twice of that for the FEM. 

For more detailed comparison between the CIP-FEM and the FEM, we consider the problem ([53]) - 
(j54p with wave number fc = 100. The traces of the CIP finite element solutions with parameters given 
by (|60p and the finite element solutions in the a;2:-plane for mesh sizes h = 1/50,1/120, and 1/200, 
and the trace of the exact solution in the xz-p\ane, are plotted in Figure [101 The shape of the CIP 
finite element solution is roughly same as that of the exact solution for h = 1/50. They match very 
well for h = 1/120 and even better for h = 1/200. While the finite element solution has a wrong 
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Figure 6: Surface plots of the real parts of the hnear mterpolant (left), the CIP finite element solution 
with parameters given by ([56| (center), and the finite element solution (right), for fc = 100 on the 
mesh with mesh size h = 1/100. 




Figure 7: h{k, 0.5) versus k for the CIP-FEM (left) with parameters given by (f56l) and for the FEM 
(right), respectively. The dotted lines give lines of slope —1.5 in the log-log scale. 

shape near the origin for h = 1/50 and h = 1/120 and only has a correct shape for h — 1/200. The 
phase error appears in all the three cases for the finite element solution. We remark that the figures 
in the left of Figure [TU] look almost the same as those in the left of Figure 6.11 in [3^, which means 
that the CIP-FEM has almost the same accuracy as the IPDG method analyzed in [3S] on the same 
mesh while using about one sixth of total DOFs of it. 

Table [1] shows the numbers of total DOFs needed for 30% relative errors in 77^-seminorm for the 
finite element interpolant, the CIP finite element solution with parameters given by ()60|) . the finite 
element solution, and the IPDG solution in [25l Subsection 6.5], respectively. The CIP-FEM needs 
less DOFs than the FEM does in all cases and much less for large wave number k. The IPDG method 
needs about six times as many total DOFs as the CIP-FEM to achieve the same accuracy but needs 
less DOFs than the FEM does when k = 100, 200, and 300. 
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Figure 8: Left graph: the relative error of the CIP finite element solution with parameters given 
by ((SD|) (solid) and the relative error of the finite element interpolant (dotted) in if^-seminorm for 
k = 5,k = 10, fc = 50, and k — 100, respectively. Right graph: the relative error of the CIP finite 
element solution with parameters given by (j60p in iJ^-seminorm computed for fc = 1, 2, • • • , 500 with 
mesh size h determined hy kh — 1 and kh = 0.5, respectively. 



k 


10 


50 


100 


200 


300 


Interpolation 


217 
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20,419 
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182,287 


CIP-FEM 
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754,507 


FEM 


397 


30,301 


229,357 


1,804,201 


6,053,461 


IPDG 


1,152 


38,088 


217,800 


1,431,432 


4,518,018 



Table 1: Numbers of total DOFs needed for 30% relative errors in i/^-seminorm for the finite element 
interpolant, the CIP finite element solution, the finite element solution, and the IPDG solution in [25l 
Subsection 6.5], respectively. 
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Figure 10: The traces of the CIP finite element solutions (left) with parameters given by (|60p and 
the finite element solutions (right) in the xz-plane for k = 100 and mesh sizes h = 1/50, 1/120, and 
1/200, respectively. The dotted lines give the trace of the exact solution in the xz-plane. 



